import numpy as np

def Ffit_fast(xb_vec, yb_vec):
    """线性拟合函数"""
    # 正向拟合
    pxy = np.polyfit(xb_vec, yb_vec, 1)
    pyx = np.polyfit(yb_vec, xb_vec, 1)
    
    # 计算正向拟合的误差
    y_fit = np.polyval(pxy, xb_vec)
    sigma1 = np.sum((yb_vec - y_fit) ** 2)
    
    # 计算反向拟合的误差
    x_fit = np.polyval(pyx, yb_vec)
    sigma2 = np.sum((xb_vec - x_fit) ** 2)
    
    if sigma2 < sigma1:
        k = 1 / pyx[0]
        b = -pyx[1] / pyx[0]
        return np.array([k, b])
    return pxy

def Fdecay_fhalf(Pt, Pf, Fs):
    """衰减计算函数"""
    Ts = 1 / Fs
    Noffset = int(np.ceil(1e-6))  # 1us
    
    # 找到RF关闭时刻
    DiffPf = np.abs(np.diff(np.abs(Pf)))
    Noff = np.argmax(DiffPf) + 1
    
    # 计算稳态振幅
    ssPtA = np.mean(np.abs(Pt)[Noff - 4*Noffset : Noff - Noffset])
    ssPtA_threshold = ssPtA * 0.5
    
    # 找到衰减结束点
    Ndecay = np.max(np.where(np.abs(Pt) > ssPtA_threshold)[0])
    
    # 计算加载Q值和fhalf
    VC = Pt
    start_idx = Noff + Noffset
    NormVCA = np.abs(VC[start_idx:Ndecay]) / np.abs(VC[start_idx])
    t2 = np.arange(len(NormVCA)) * Ts
    
    p = Ffit_fast(t2, np.log(NormVCA))
    return -p[0]